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1.  Introduction  to  Post-Stall  Aerodynamics 

Fighter  aircraft  have  evolved  for  decades  toward  high  speeds  and 
high  wing  loadings.  Advances  in  air-to-air  weapons  are  bringing  that 
line  of  evolution  to  a  halt,  and  the  advantage  is  passing  to  fighter 
with  lower  wing  loadings  and  great  agility.  The  new  design  direction  is 
apparent  in  the  F-16XL  shown  in  figure  1.1.  The  Rockwell  model  in 
figure  1.2  displays  further  concessions  to  maneuverability. 

Herbst  (1980)  has  introduced  the  concept  of  "supermaneuverability" 
to  characterize  the  design  of  future  fighters.  Supermaneuverability  has 
two  aspects:  decoupling  of  orientation  from  flight  path  to  allow  the 
aircraft  to  function  as  a  rotating  gun  platform,  and  flight  at  angles  of 
attack  far  beyond  the  limits  of  stall. 

Herbst  uses  simulated  fighter  engagements  to  assess  the  merits  of 
post-stall  (PST)  maneuvers.  Figure  1.3  shows  time  histories  of  speed  V 
and  angle  of  attack  a  for  a  PST  fighter  engaging  a  normal  adversary. 
The  engagement  begins  with  the  fighters  approaching  each  other  at  a 
closing  speed  of  200  m/sec.  The  PST  fighter  executes  a  180°  turn  at  an 
of  90*3,  broadside  to  his  flight  path.  His  engine  thrust  balances 
centrifugal  force  as  his  flight  speed  drops  to  20  m/sec.  Realigning  his 
thrust  with  his  new  flight  path,  he  accelerates  up  the  flight  path  of 
his  hapless  adversary. 

Flight  at  such  large  angles  introduces  phenomena  outside  the 
boundaries  of  traditional  aerodynamics.  Some  of  those  phenomena  are 
evident  in  figure  1.4: 

•  Flow  separates  all  around  the  planform  edges.  This  produces  the 
bag-shaped  vortex  sheet  shown  in  cross  section  in  figure  1.4. 
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•  Viscous  effects  are  unimportant;  a  consequence  of  the  violent 
and  unambiguous  separation  on  all  salient  edges  of  the  planform. 

•  Pressure  force  is  normal  to  the  planform  owing  to  a  complete 
loss  of  leading-edge  suction. 

•  The  vortex  sheet  shed  from  the  aircraft  is  strong,  wi  th  a 
velocity  jump  comparable  to  the  speed  V  of  the  aircraft. 

•  Shed  vorticity  interacts  strongly  with  the  aircraft,  probably 
creating  "vortex  lift"  like  vorticity  shed  from  leading  edges  of 
del ta  wings. 

•  The  vorticity  around  the  aircraft  is  nonsteady.  Kelvin- 
Helmholtz  instability  of  the  vortex  sheet  precludes  *he 
attainment  of  a  steady  Helmholtz  (free-streamline)  flow,  with 
dead  air  inside. 

The  phenomena  above  are  beyond  the  reach  of  aerodynamic  theories 
predicated  upon  quasi -steady  flow  or  small  perturbations.  Numerical 
methods,  as  described  here,  almost  certainly  are  required  to  handle  the 
nonsteady,  nonlinear  vortex  flow  surrounding  a  PST  aircraft  and  to 
predict  the  resulting  aerodynamic  forces. 
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2.  Theoretical  Foundations 

a.  Theorems 

Our  study  of  post-stall  aerodynamics  begins  with  the  Euler 
equations  for  inviscid  incompressible  flow: 

.3u  V 

p  —  +  il  •  Vu)  +  VP  =  f  .  (2.1) 

'at  ' 

V  •  u  =  0  .  (2.2) 

[n  these  equations  p  is  density,  p  is  pressure,  lu  is  velocity,  and 
f  is  a  force  per  unit  volume  exerted  on  the  fluid  by  some  external 
agency.  The  continuity  equation  (2.2)  presumes  that  p  is  constant, 
an  adequate  approximation  for  low-speed  maneuvers  though  not  for 
transonic  flight.  The  assumption  of  inviscid  flow  is  warranted  only  if 
we  can  devise  side  conditions  that  accurately  represent  the  conversion 
of  boundary-layer  vorticity  into  free  vorticity  at  planform  edges,  a 
matter  discussed  in  part  (c)  of  this  chapter. 

Equations  (2.1)  and  (2.2)  apply  to  fluid  outside  the  aircraft 
surface,  on  which  normal  velocity  boundary  conditions  are  imposed. 
Conceptually  and  computationally,  it  is  often  convenient  to  regard  the 
fluid  as  filling  all  space,  with  (2.1)  and  (2.2)  applying  everywhere, 
and  with  the  aircraft  surface  represented  by  some  distribution  of 
singularities,  sources  and  bound  vortex  sheets  being  popular  choices. 
Completely  equivalent  to  a  vortex  sheet  is  a  dipole  sheet,  which  is  a 
volumetric  force  concentrated  on  the  surface  of  the  body  and  distributed 
over  the  surface  so  as  to  enforce  the  normal  velocity  boundary 
condition.  Such  a  representation  is  achieved  by  writing  £  as 
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where  Ap  is  the  pressure  jump  across  the  dipole  sheet,  6(n)  is  a 
delta  function  of  the  normal  coordinate  n,  and  n  is  the  unit  outward 
normal  vector.  The  pressure  jump  Ap  is  a  function  of  time  and  of  two 
surface  coordinates  and  is  tailored  to  insure  that 

u  •  n  =  Vp  .  (2.4) 

where  vp  is  the  local  normal  component  of  boundary  velocity.  With 
this  representation,  equations  (2.1)  and  (2.2)  apply  throughout  a 
hypothetical  unbounded  fluid  pervading  the  aircraft  body. 

Recall  that  ^  is  a  volumetric  force  exerted  ^  an  external  agency 
(the  aircraft)  o^  the  fluid,  the  net  force  exerted  by  the  fluid  on  the 
aircraft  is  therefore  given  by  the  volume  integral 


V 


where  the  region  of  integration  V  extends  over  all  space.  A  similar 
formula  provides  the  net  moment  exerted  by  the  fluid  on  the  aircraft: 


X  X  f  dV  .  (2.6) 

V 

The  dipole  representation  (2.3)  allows  an  immediate  evaluation  of  £ 
and  M  without  a  need  for  computing  pressure. 

There  is  no  simple  integral  relation  between  F  or  M  and  u, 
because  £  does  not  fall  toward  zero  fast  enough  to  permit  a  proper 
evaluation  of  integrals  like 
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There  are,  however,  simple  relations  of  £  and  M  to  vorticity  w,  the 
curl  of  the  velocity  field; 

03  =  V  X  u  .  (2.7) 

Vorticity  U3  is  a  complete  flow  representation  when  £  is  divergence 
free.  Velocity  may  then  be  written  as  the  curl  of  a  vector  potential  , 


u  »  7  X  1(3 


(2-.8) 


with 


“  “  0)  •  (2.9) 

Given  u,  one  solves  the  Poisson  equation  (2.9)  for  from  which  u 
follows  by  (2.8).  Of  course  w  includes  both  actual  "free"  vorticity 
o3p  in  the  fluid,  and  "bound"  vorticity  wj,  associated  with  the  dipole 
sheet  representing  the  body: 


£  =  ub  +  wf  .  (2.10) 

£b  is  the  concentrated  vortex  sheet  that  would  appear  in  a  vortex-p?nel 
code. 


Vorticity  evolves  according  to  an  equation  obtained  by  taking  the 
curl  of  (2.1); 


903 

--  +  U  •  703  -  W  •  Vu  =  7  X  g  (2.11) 

3t  ~ 
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where 


g  =  f/p  .  (2.12) 

We  can  derive  a  relation  between  F  and  x  x  to  by  taking  the  cross- 
product  of  £  with  (2.11)  and  integrating  over  all  space,  using  the 
facts  that  ui  ^nd  g  are  zero  outside  some  finite  region.  The  net 
moment  M  bears  a  corresponding  relation  witn  x  x  x  x  co.  Appendix  A 

rw 

shows  how  the  integral  relationships  are  derived. 

The  force  and  moment  theorems  must  be  derived  separately  for  two 
and  three-dimensional  flows,  because  a  two-dimensional  vorticity  field 
is  unbounded  in  the  third  dimension,  which  changes  the  convergence 
properties  of  integrals.  A  two-dimensional  field  of  vorticity  aligned 
with  the  2-axis  has  tne  form 

01  *  u>(x.  y,  t)  Iz  .  (2.13) 

Net  force  per  unit  length  in  the  z-di recti  on  has  the  form 

d 

F  =  -  p- 

d- 

A 

where  the  integration  extends  over  the  full  (x,  y)  plane,  or  at  least 
over  every  part  containing  non-zero  vorticity.  Moment  per  unit  length 
has  the  form 


X  (uer)  dA 


(2.14) 


M 


0  d 
2  dt 


y*  {x2  +  y2)  (ue^)  dA 
A 


(2.15) 


directed  along  the  z-axis.  Those  are  reactions  exerted  on  the  aircraft 
by  the  fluid. 
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The  corresponding  results  in  three  dimensions  are  as  follows: 


F 


P 

2  dt 


I 


X  X  u)  dV 


(2.16) 


M 


P  d 
3  dt 


(£  X  U  )  dV 


(2.17) 


Equations  (2.14),  (2.15)  or  (2.16),  (2.17)  allow  us  to  evaluate  forces 
and  moments  from  vorticity  without  the  intervention  of  pressure. 
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b.  Applications  in  Two  Dimensions 

Figure  2.1  shows  a  typical  situation  where  we  might  want  to  apply 
the  theorems  of  part  (a)  in  two  dimensions.  A  thin  airfoil  starts  from 
rest,  shedding  vorticity  from  its  trailing  edge  and  possibly  from  its 
leading  edge.  The  airfoil  can  be  represented  as  a  sheet  of  bound 
vortici ty , 


y,  t)  *  Y(s,  t)  6(n)  ,  (2.18) 

where  s  and  n  are  coordinates  parallel  and  normal  to  the  airfoil 
respectively,  5(n)  is  the  delta  function  expressing  the  concentration 
of  the  bound  vortex  sheet,  and  Y(s,  t)  is  the  velocity  jump  from  the 
bottom  to  the  top  of  the  airfoil: 


Y  =  U2  -  ui 


(2.19) 


According  to  the  Kelvin  circulation  theorem,  the  net  vorticity  over 
the  entire  plane  remains  zero  for  all  time  if  the  flow  starts  from  rest: 


r 


dA 


0 


(2.20) 


Equation  (2.20)  is  an  important  constraint  on  the  assignment  of  bound 
vorticity  to  represent  the  airfoil.  The  drag  and  lift  on  the  airfoil 
follow  fr  n  (2.14): 


D 


dt 


/ 


-  p —  /  y(u)b  +  ojf)  dA  , 


J'  x(a)b  +  Uf)  dA 

A 


(2.21) 


(2.22) 


d 

L  =  0  — 
dt 
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Figure  2,1.  Bound  and  free  vorticity  in  two  dimensions 
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If  the  airfoil  eventually  stops  shedding  vorticity,  then  the  free 
vorticity  tof  wraps  into  a  starting  vortex  with  constant  circulation 
and  a  centroid  fixed  in  space,  i.e. 


d 

dt 


/ 


>c  ujf  dA  -►O 


(2.23) 


The  bound  vorticity  becomes  steady  in  coordinates  (X  =  x  +  Vt,  y) 
moving  leftward  with  the  airfoil  at  its  velocity  V.  Thus 


D 


P—  yu>b(X,  y)  dXdy  =  0, 


(2.24) 


and 


L 


pV 


If 

f 


(X  -  Vt)  ojbfx,  -y)  dXdy 


ojb  dA  =  pVTb  , 


A 


(2.25) 


where  Tb  is  the  circulation  around  the  airfoil  defined  in  the  clock¬ 
wise  sense  familiar  in  aerodynamics.  Equation  (2.15)  is  D'Alembert's 
paradox  of  zero  drag  in  steady  flow,  and  (2.16)  is  the  Kutta-Joukowski 
theorem  for  lift.  Those  results  should  help  build  confidence  in  (2.14) 
as  the  general  force  law  for  non-steady  two-dimensional  flow. 


A  key  assumption  in  our  derivation  of  the  Kutta-Joukowski  theorem 
was  that  vortex  shedding  ceases  after  some  time,  so  the  flow  becomes 
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steady  in  airfoil -fixed  coordinates.  Violation  of  that  assumption  is 
one  of  the  defining  characteristics  of  post-stall  aerodynamics. 

One  way  vortex  shedding  can  continue  indefinitely  is  for  shedding 
to  occur  at  both  the  leading  and  trailing  edges  of  the  airfoil.  As  long 
as  vorticity  shed  from  the  leading  edge  has  sign  opposite  that  shed  from 
the  trailing  edge,  shedding  can  continue  without  violating  the  vorticity 
conservation  theorem  (2.20).  As  the  two  vortex  sheets  lengthen,  they 
might  even  assume  an  appearance  of  steadiness  in  the  vicinity  of  the 
ai rfoi 1 . 

Helmholtz  developed  the  holegraph  method  to  study  such  flows  around 
flat  plates,  as  illustrated  in  figure  2.2(a)  (Lamb  1932).  Helmholtz 
derived  the  formula 

ir  sin  a 

F  =  -  pV2c  (2.26) 

4  +  IT  sin  a 

for  force  per  unit  length  normal  to  an  infinite  flat  plate  (lamina)  of 
chord  c.  We  might  have  hoped  that  Helmholtz’s  solutions  become 
asymptotically  valid  for  post-stall  flows  at  long  times  after  pitch-up. 
Such  hope  is  unwarranted,  because  the  free  shear  layers  of  Helmholtz 
flows  are  subject  to  a  Kel vin-Helmholtz  instability  as  shown  in  figure 
2.2(b).  Steady-state  free-streamline  theory  predicts  a  zero-velocity 
wake  with  a  uniform  pressure  equal  to  the  ambient  pressure  far  upstream. 
Eddies  resulting  from  the  instability,  however,  scavenge  the  wake  and 
lower  the  backpressure  on  the  plate. 

Equation  (2.26)  predicts  a  drag  coefficient  of  0.88  for  a  plate 
normal  to  the  flow,  whereas  the  measured  value  is  1.98  (Hoerner  1965). 
The  125  percent  increase  is  attributable  to  nonsteady  vortex  suction, 
not  too  surprising  in  view  of  the  origin  of  aerodynamic  forces  evident 
in  (2.21)  and  (2.22). 
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(o)  HELMHOLTZ  FLOW 


Figure  2.2.  Helmholtz  flow  (a)  and  its  nonsteady  counterpart 
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Vortex  suction  plays  an  important  role  in  the  lift  of  delta  wings 
at  large  angles  of  attack.  Indeed  the  flow  around  a  delta  wing  is  a 
prototype  for  post-stall  aerodynamics.  Figure  2.3  shows  a  delta  wing 
inclined  at  a  substantial  angle  a  to  the  flow.  A  plane  A  traverse 
to  the  delta-wing  planform  moves  with  the  fluid  at  speed  V.  The  delta 
wing  intersects  the  transverse  plane  as  a  slit  of  varying  span  b(t), 
moving  downward  at  a  speed  V  sin  a.  Vorticity  shed  from  the  leading 
edge  of  the  delta  wing  appears  as  free  vorticity  in  the  transverse 
plane.  If  the  delta  wing  is  slender,  then  the  flow  in  the  transverse 
plane  can  be  treated  like  the  two-dimensional  time-dependent  flow  of 
figure  2.1.  In  particular,  the  lift  per  unit  length  along  the  delta 
wing  is  given  by  equation  (2.22)  evaluated  in  the  fluid-fixed  transverse 
plane. 

There  is  no  fully  satisfactory  theory  for  the  vortex  lift  of  delta 
wings.  Vortex-lift  calculations  rely  instead  on  analogies.  Hoerner 
(1975)  describes  an  end-plate  analogy,  whereby  the  leading-edge  vortices 
are  regarded  as  acting  like  end  plates,  increasing  the  effective  span  in 
the  transverse  plane.  Polhamus  (1971)  achieved  good  correlation  with 
data  on  the  basis  of  a  leading-edge  suction  analogy,  where  the  leading- 
edge  suction  that  reduces  the  drag  of  an  attached  flow  is  rotated  90°  to 
produce  lift  normal  to  the  planform.  Neither  analogy  is  convincing  as  a 
basis  for  a  venture  into  post-stall  aerodynamics. 

A  key  difficulty  with  the  concept  of  vortex  lift  is  evident  from 
equation  (2.22):  vortex  lift  arises  from  the  time  dependence  of  free 
vorticity  in  the  two-dimensional  flow  of  figure  2.1  or  in  the  transverse 
plane  of  figure  2.3.  Steady  vortices  in  such  situations  exert  no  net 
force  on  bodies.  A  satisfactory  method  for  predicting  the  vortex  lift 
of  delta  wings  should  be  among  the  first  products  of  post-stall 
aerodynamics. 
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c.  Vortex  Shedding 

We  have  assumed  so  far  that  the  flow  is  inviscid,  though  it  may 
contain  "free”  vorticity  and  "bound"  vorticity  representing  the  surfaces 
of  bodies.  The  distinction  between  bound  and  free  vorticity  is  peculiar 
to  the  assumption  of  inviscid  flow  and  would  disappear  if  we  were 
willing  to  deal  with  viscous  diffusion.  "Bound"  vorticity  is  simply 
vorticity  in  boundary  layers,  subject  to  viscous  and  inertial  dynamics 
about  equally.  "Free"  vorticity  is  vorticity  that  has  escaped  far 
enough  from  boundaries  to  be  subject  mainly  to  inertial  rather  than  to 
viscous  dynamics.  The  escape  process  is  called  vortex  shedding,  meaning 
the  transformation  of  bound  into  free  vorticity. 

The  reason  for  maintaining  a  distinction  between  bound  and  free 
vorticity  is  that  we  hope  to  analyze  post-stall  aerodynamics  using 
inviscid  flow  equations.  This  hope,  however,  requires  us  to  postulate  a 
rule  for  transforming  bound  into  free  vorticity  without  explicit 
reliance  on  the  Navier-Stokes  equations  of  viscous  flow.  Such  a  rule  is 
unlikely  to  become  available  for  bodies  of  general  shape,  where  both  the 
location  and  rate  of  vortex  shedding  are  variable.  We  therefore  confine 
attention  to  sal ient  edges ,  from  which  vorticity  can  be  expected  to  shed 
almost  all  the  time. 

The  search  for  inviscid  shedding  rules  is  a  challenge  for 
computational  fluid  dynamics.  Some  authors  (e.g.  Hitzel  and  Schmidt 
1984)  invoke  compressibility  as  a  shedding  mechanism.*  Compressibility 
will  be  of  no  help  to  us,  because  we  want  to  study  incompressible  flows. 
It  is  far  from  clear,  moreover,  that  compressible  effects  can  mimic 
viscous  flow  separation. 


*As  did  Rayleigh,  interestingly,  shortly  before  Prandtl's  1905  paper  on 
the  boundary-layer  concept. 
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Another  school  of  thought  holds  that  finite  difference  schemes  by 
their  very  nature  adequately  simulate  the  effects  of  viscosity  with 
respect  to  vortex  shedding  at  salient  edges  (Rizzi  and  Eriksson  1984). 
There  is  merit  to  this  viewpoint  if  velocity  is  one  of  the  advected 
quantities.  Upwind  differencing,  in  particular,  seems  almost  certain  to 
project  velocity  discontinuities  (vortex  sheets)  from  a  sharp  edge  into 
the  computational  grid.  We  have  based  the  two-dimensional  computations 
presented  here  on  the  vorticity  equation,  however,  which  sustains  zero 
vorticity  as  a  solution  unless  free  vorticity  is  injected  deliberately 
into  the  finite-difference  grid.  We  must  provide  an  explicit  algorithm 
for  the  the  creation  of  free  vorticity. 

Our  vortex  shedding  algorithm  is  based  on  the  fact  mentioned  above 
that  bound  vorticity  is  vorticity  in  boundary  layers.  A  boundary  layer 
cannot  flow  around  a  sharp  edge  but  must  separate  from  it  if  the 
direction  of  flow  is  from  the  material  boundary  toward  free  fluid.  The 
rate  of  creation  of  free  vorticity  is  the  integrated  flux  of  boundary- 
layer  vorticity  flowing  across  the  edge. 

Figure  2.4(a)  illustrates  the  principles  of  vortex  shedding  in  two 
dimensions  when  the  flow  is  directed  away  from  the  body  both  above  and 
below  an  edge.  Consider  the  situation  above  the  edge,  and  adopt  local 
coordinates  (x,  y)  such  that  x  is  parallel  to  the  upper  surface  and 
y  is  normal  to  it.  Vorticity  can  be  written  everywhere  as 

0)  =  9V/3X  -  3u/3y  .  (2.27) 

In  the  boundary  layer,  v  is  small  and  x-deri vati ves  are  small,  so  on 
both  counts  vorticity  can  be  approximated  as 

u)  »  -  3u/3y  .  (2.28) 
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The  net  vorticity  flux  across  the  edge  in  the  upper  boundary  layer  is 


5 


(u-Uq)  03  dy 


(2.29) 


where  6  is  the  depth  of  the  boundary  layer,  and  Uq  accounts  for  the 
possibility  that  the  upper  surface  of  the  body  is  moving  at  a  speed  Ug 
parallel  to  itself.  Adopting  the  boundary-layer  approximation  (2.28) 
for  03  and  integrating  through  the  boundary  layer,  we  obtain 


2 

as  the  expression  for  net  vorticity  flux  consistent  with  thin  boundary 
layers.  Here  u^  is  the  value  of  u  at  the  top  of  the  boundary  layer, 
i  .e.  the  surface  slip  velocity  of  inviscid  flow  theory,  while  Ug  is 
the  value  of  u  at  the  boundary  consistent  with  the  actual  no-slip 
boundary  condition. 


/3u 

(U-Ug)  — 

3y 


dy  =  - 


Though  (2.30)  is  based  on  boundary-layer  theory,  the  formula 
contains  no  explicit  reference  to  viscosity.  This  supports  the  hope 
that  post-stall  aerodynamics  can  be  treated  by  inviscid  flow  methods  as 
long  as  the  shedding  rule  correctly  embodies  the  boundary-layer  concept. 

Similar  reasoning  about  the  underside  of  the  surface  leads  to  the 
total  vorticity  flux  expressed  in  figure  2.4(a).  Notice  that 


(U2  -  Uo)2  (ui  -  Uo)2 
2  2 


(U2  -  ui) 


(2.31) 


*  (ii  -  UglY  , 
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Figure  2.4.  Transformation  of  bound  into  free  vorticity. 
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i  .e.  the  net  vorticity  flux  equals  the  strength  y  of  the  vortex  sheet 
at  the  edge  times  the  free-stream  velocity  relative  to  the  edge  averaged 
above  and  below  the  lifting  surface. 

The  situation  in  figure  2.4(a)  is  not  the  only  possibility.  The 
free-stream  velocity  relative  to  the  edge  might  be  directed  inward  from 
the  edge  above  the  surface,  below  the  surface,  or  both.  Figure  2.4(b) 
illustrates  the  consequences  of  an  inward-directed  velocity  above  the 
surface.  No  boundary-layer  vorticity  sheds  from  the  upper  surface,  and 
only  the  first  term  on  the  left  of  (2.31)  contributes  to  a  vorticity 
flux.  The  simple  interpretation  of  vorticity  transformation  afforded  by 
the  right  of  (2.31)  is  lost,  but  the  phenomenology  is  clear. 

Our  vortex  shedding  rule  must  specify  not  only  how  much  free 
vorticity  is  created,  but  also  where  it  is  place  in  the  computational 
grid.  Figure  2.5  shows  how  the  placement  is  effected  in  the  code 
described  in  later  sections.  Free  vorticity  is  injected  into  the 
nearest  grid  cell  lying  entirely  downwind  of  the  surface  from  which 
vorticity  is  being  shed. 

There  is,  of  couse,  a  certain  arbitrariness  in  this  prescription. 
We  trust  that  the  vortex  shedding  process  is  self-equilibrating,  such 
that  shed  vorticity  erases  singularities  near  edges  and  proceeds 
independently  of  fine  details  near  the  edges.  Viscous  flow  analyses  or 
comparison  with  experiments  will  be  needed  to  validate  that  trust. 
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3.  Numerical  Methods  for  Two-dimensional,  Unsteady,  Separated  Flow 

a.  Introduction 

Methods  for  numerical  modeling  of  lift-generated  vortex  flows  can 
be  roughly  divided  into  those  methods  that  employ  discrete  elements  for 
free  vorticity,  as,  for  example,  Geising  (1968),  and  Hoeijmaker  and 
Vaatstra  (1983)  and  those  that  represent  free  vorticity  as  a  continuous 
field,  as,  for  example,  Rizzi  and  Eriksson  (1984).  We  have  chosen  to 
represent  vorticity  as  a  continuous  field  largely  because  we  had  much 
greater  confidence  in  the  stability,  accuracy,  and  physical 
reasonableness  of  solutions  produced  by  Euler  solvers  than  in  those 
produced  by  discrete  element  methods. 

Given  the  choice  of  a  continuous  field  representation  for  free 
vorticity,  a  complete  predictive  code  for  two-dimensional,  unsteady, 
separated  flow  will  have  three  ingredients: 

1.  A  means  for  representing  the  effect  of  a  solid  obstacle  on  the 
flow  velocity, 

2.  A  means  for  shedding  vorticity  from  the  body  into  the  free 
vorticity  field,  and 

3.  A  means  for  advecting  the  free  vorticity. 

The  third  ingredient  can  be,  and  for  our  work  is,  a  garden-variety  Euler 
solver.  The  first  and  third  ingredients  may  be  combined  in  a  code  that 
employs  body-fitted  coordinates,  but,  as  we  shall  see,  it  is  possible  to 
build  a  code  that  models  arbitrary  bodies  without  the  complexities  of 
body  fitted  coordinates.  The  third  ingredient  being  the  foundation  upon 
which  the  whole  is  built,  we  proceed  with  a  brief  description  of  the 
Euler  solver. 
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b.  Advection  of  Free  Vorticity 

We  have  taken,  as  the  starting  point  for  our  work,  the  Poseidon 
Hydro  code,  developed  originally  in  1976  for  the  simulation  of  two 
dimensional  vortex  flows  in  stratified  fluid.  The  code  has  been  exten¬ 
sively  exercised,  validated,  and  refined  under  a  variety  of  programs  for 
DARPA  and  for  the  U.  S.  Navy,  but  the  only  unclassified  description  of 
the  code  and  its  applications  is  given  by  St.  Cyr  (1980).  The  code  has, 
over  the  years,  acquired  many  special  features  directed  toward  the 
programs  for  which  it  was  exercised.  Our  first  step  under  the  current 
program  was  to  gather  only  the  essential  nucleus  of  this  software  into 
an  entirely  new  code.  We  used  an  established  test  case  (a  vortex  pair 
in  unstratified,  inviscid,  incompressible  fluid)  to  verify  that  the 
"new"  code  reproduced  exactly  the  unstratified,  inviscid  dynamics  of  the 
Hydro  code.  This  done,  we  had  a  code  which  had  all  the  advantages  of 
many  years  ot  testing  behind  it,  but  without  the  breadboard  additions 
that  old  codes  inevitably  acquire. 

The  Hydro  code,  and  its  descendant,  LIFTER,  are  formulated  in  terms 
of  vorticity  and  streamfuncti on .  This  formulation  has  many  nice 
features  for  two  dimensional  flow  and  is  by  far  the  most  popular  formu¬ 
lation  for  two-dimensional  incompressible  flow  codes.  A  code  formulated 
in  terms  of  vorticity  and  streamfuncti  on  has  just  two  essential  ingre¬ 
dients:  a  scheme  for  transporting  vorticity,  and  a  scheme  for  computing 
the  streamfuncti on  given  the  vorticity. 

The  advective  transport  scheme  is  a  derivative  of  the  original 
SHASTA  scheme  of  Boris  and  Book  (1973).  SHASTA  is  monotonic  for  time 
steps  such  that  the  local  Courant  number  is  always  less  than  0.5.  That 
is,  transport  never  produces  new  local  extrema  in  transported  field 
quantities.  Monotonicity  is,  in  fact,  the  key  design  requirement  for  a 
whole  class  of  schemes  which  have  come  to  be  called  flux-corrected 
transport  schemes,  of  which  SHASTA  is  an  example.  SHASTA  is  as  close  to 
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being  second  order  accurate  in  space  as  it  is  possible  to  be  without 
sacrificing  monotonicity  and  exhibits  exceptionally  small  residual 
numerical  diffusion.  Most  schemes  that  are  formally  second  order 
accurate  and  stable  are  substantially  more  diffusive  than  SHASTA.  The 
algorithm  that  is  implemented  in  LIFTER  differs  from  the  original  SHASTA 
only  in  that  a  slightly  more  accurate  flux  correction  term  is  employed 
in  LIFTER.  Time  differencing  is  done  using  a  second  order  Runge-Kutta 
(predictor-corrector)  scheme. 

The  streamfunction  is  computed  from  vorticity  by  solving  a  Poisson 
equation : 


3  _  (jJ 

Two  issues  must  be  dealt  with  in  solving  this  equation:  representation 
of  an  infinite  domain  on  a  finite  grid  and  an  efficient  way  of  solving 
the  difference  equations.  If  all  of  the  voriticity  is  captured  within 
the  computational  domain,  an  infinite  domain  may  be  represented  by 
computing  the  value  of  the  streamfunction  on  the  boundary  of  the 
computational  domain  using  a  discrete  form  of  the  Biot-Savart  integral 
for  the  streamfunction,  and  then  solving  the  Poisson  equation  subject  to 
the  computed  Oirichlet  boundary  condition.  This  reasonable  strategy  was 
employed  in  the  Hydro  Code  as  described  by  St.  Cyr  (1980).  We  have 
subsequently  adopted  the  approach  of  solving  the  Poisson  Problem  on  a 
very  large  grid  using  fast  Fourier  transforms.  Despite  the  fact  that 
the  computational  domain  must  be  substantially  enlarged  for  this  latter 
strategy,  it  is  significantly  faster  than  the  original  .  The  speed 
advantage  arises  because  the  operations  count  for  the  Biot-Savart 
integral  is  0(n3)  while  the  operations  count  for  the  Fast  Fourier 
Transform  is  0(N  In  N).  Furthermore,  the  Block  Cyclic  Reduction 
algorithm  formerly  used,  while  having  and  operations  count  that  is  0(N 
In  N),  involves  tri diagonal  sweeps  that  do  not  vectorize  on  the  CRAY. 
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c.  Representation  of  Wing  Sections  by  Bound  Vorticity 

The  velocity  field  computed  by  the  free  vorticity  code  is  different 
from  the  actual  velocity  field  around  a  solid  body  by  a  flow  which  is 
everywhere  irrotational  except  within  the  solid  body.  For  our  purposes, 
the  solid  body  will  always  be  a  one-dimensional  curve  within  the  r^mpu- 
tational  domain.  In  order  to  compute  the  actual  velocity  field,  we  must 
obtain  a  flow  that  is  irrotational  everywhere  but  along  the  curve  and 
that  cancels  along  the  curve  the  normal  velocity  due  to  free  vorticity. 

There  are  a  variety  of  techniques  available  to  compute  the  required 
velocity  field,  but  by  far  the  most  attractive  are  the  so-called  panel 
methods.  These  methods  compute  the  required  velocity  field  as  if  it 
resulted  from  a  distribution  of  flow  field  singularities  distributed 
along  the  body.  The  great  power  of  these  methods  derives  from  the  fact 
that  the  body  is  represented  by  polygonal  elements,  permitting  the 
representation  of  nearly  arbitrary  geometries. 

Of  the  many  possible  choices  for  the  flow  field  singularities,  by 
far  the  most  natural  for  the  problem  at  hand  is  a  sheet  of  vorticity. 
Vortex  sheet  panel  codes,  as  in  Raj  and  Gray  (1978),  are  a  standard 
technique  for  the  aerodynamics  of  thin  airfoils. 

Figure  3-1  illustrates  most  of  the  essential  features  of  the  panel 
code  that  has  been  embedded  into  lifer.  The  wing  is  represented  as 
series  of  N  panels  of  arbitrary  length  and  orientation,  each  panel  being 
a  sheet  of  vorticity  of  constant,  unknown  strength.  Panel  a  induces  a 
normal  velocity  at  a  chosen  control  point  on  panel  S;  the  value  of  this 
velocity  for  unit  vorticity  per  unit  length  on  panel  a  is  the 
influence  coefficient  Cag.  The  sum  of  all  of  the  induced  velocities  at 
panel  6  must  exactly  cancel  the  normal  velocity  at  panel  S  due  to 
the  free  vorticity.  We  thus  have  N  equations  for  the  N  unknown 
panel  vortex  strengths. 
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The  matrix  of  influence  coefficients  computed  by  the  preceding 
prescription  will  be  numerically  ill-conditioned,  as  a  consequence  of 
being  nearly  singular.  This  near  singularity  is  a  consequence  of  the 
exact  singularity  of  the  continuum  integral  equation  being  approximated. 
The  singularity  of  the  integral  equation  is  reflection  of  the  well-known 
fact  that  the  flow  around  an  airfoil  can  have  an  arbitrary  net 
circulation.  This  indeterminacy  is  resolved  in  conventional 
aerodynamics  by  insisting  that  a  stagnation  point  coincide  with  the 
sharp  trailing  edge  of  an  airfoil.  This  requirement  is  commonly  known 
as  the  Kutta  condition. 

If  the  flow  around  the  airfoil  is  permitted  to  separate,  the 
arguments  that  lead  to  the  Kutta  condition  fall  apart.  From  the 
absolutely  fundamental  Kelvin  circulation  theorem,  we  know  that  the 
total  circulation  of  the  flow  must  be  preserved.  If  we  start  with  an 
airfoil  at  rest,  then  the  indeterminacy  is  resolved  by  insisting  the 
total  circulation  of  the  flow  be  zero.  This  condition,  which  we  have 
chosen  to  call  the  Kelvin  condition,  is  used  to  replace  one  of  the  rows 
in  the  matrix  of  influence  coefficients. 

All  that  has  been  said  so  far  is  reasonable  and  even  obvious,  but 
some  details  have  been  glossed  over.  For  example,  the  most  obvious 
choice  of  control  points,  the  center  of  each  panel,  produces  an 
oscillating  distribution  of  bound  vorticity  on  a  flat  plate  in  a  uniform 
free  stream.  Moving  the  control  point  to  the  quarter  chord  of  each 
panel  results  in  a  vorticity  distribution  that  more  nearly  approximates 
the  know  solution  to  the  continuum  equations.  This  naturally  begs  the 
question  as  to  why  the  quarter  chord  and  not  the  three  quarter  chord? 
Or,  to  put  it  another  way,  what  is  the  difference  between  the  quarter 
chord  and  the  three  quarter  chord  for  a  flat  plate  at  90  degrees  angle 
of  attack?  Similarly,  which  of  the  equations  is  to  be  replaced  by  the 
Kelvin  condition? 
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We  have  alleviated  the  arbitrariness  of  the  choice,  at  least  to 
some  extent,  by  solving  for  the  bound  vorticity  with  one  arbitrary 
choice  of  control  points  and  for  its  mirror  image,  and  then  averaging 
the  results.  The  Kelvin  condition  in  each  case  replaces  a  control  point 
equation  at  one  end  of  the  panel.  While  the  choices  are  still 
arbitrary,  the  procedure  at  least  preserves  any  symmetries  the  problem 
may  have  had  before  descretization. 

A  further  detail  we  must  deal  with  is  that  once  we  have  the 
singularity  distribution,  we  must  compute  the  velocity  field  due  to  the 
singularity  distribution  at  every  point  in  the  computational  grid. 
Calculating  this  velocity  field  from  a  Biot-Savart  integral  is  imprac¬ 
tical,  not  only  because  of  the  operation  count,  but  also  because  an 
impracticably  large  number  of  influence  coefficients  would  have  to  be 
stored.  The  most  economical  way  to  generate  a  velocity  field  from  the 
bound  vorticity  is  to  assign  the  bound  vorticity  to  grid  points  and 
solve  a  Poisson  equation  as  for  the  streamfunction  due  to  free  vorti¬ 
city.  This  seemingly  crude  procedure  seems  to  yield  remarkably  good 
results  as  will  be  evident  from  the  streamfunction  plots  in  Chapter  4. 

d.  Conversion  of  Bound  Vorticity  into  Free  Vorticity 

In  Chapter  Two,  we  discussed  at  length  the  physics  of  our  model  for 
vorticity  shedding.  Here  it  remains  only  to  report  on  the  details  of 
the  method  by  which  this  shedding  model  is  implemented.  The  mean  velo¬ 
city  at  the  wingtip  is  computed  by  interpolating  the  velocity  field  on 
the  grid,  and  subtracting  the  velocity  of  the  wingtip,  if  the  wingtip  is 
not  stationary  in  the  grid.  The  vorticity  on  the  end  panel  is  presumed 
not  to  contribute  to  the  mean  velocity  at  the  wingtip.  The  mean  tangent 
velocity  is  computed  by  dotting  the  mean  velocity  with  an  outward  facing 
unit  tangent  vector.  The  tangent  velocity  on  one  side  of  the  wing  is 
taken  to  be  the  mean  tangent  velocity  plus  y/2  and  the  tangent 
velocity  on  the  other  side  of  the  wing  is  taken  to  be  the  mean  tangent 
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velocity  minus  Y/2,  where  Y  is  the  bound  vorticity  at  the  wingtip. 
Given  the  velocity  on  either  side  of  the  wing,  it  is  a  straightforward 
matter  to  compute  the  rate  of  shedding  from  the  formulas  in  figure  2.4. 

Once  the  rate  of  shedding  has  been  computed,  the  vorticity  is  added 
to  the  free  vorticity  in  the  grid.  In  order  to  ensure  that  the  newly 
freed  vorticity  is  not  trapped  at  the  wingtip,  the  vorticity  is  added  to 
a  control  volume  that  does  not  contain  the  wingtip. 

This  description  makes  clear  a  possible  weak  point  in  the  model. 
The  rate  of  vorticity  shedding  is  proportional  to  the  value  of  the 
vorticity  on  the  last  panel.  For  an  impulsively  started  flat  plate,  the 
vorticity  at  the  end  of  ihe  plate  is  infinite,  and  yet  the  rate  of 
vorticity  shedding  computed  by  our  model  is  dependent  upon  the 
artificially  finite  value  computed  as  a  result  of  the  discretization 
process.  The  only  reason  to  believe  that  such  an  arbitrary  procedure 
could  yield  correct  results  is  that  the  process  is  self-correcting. 
That  is,  one  effect  of  shedding  vorticity  is  to  relieve  the  singularity 
at  the  wingtip.  If  too  small  an  amount  of  vorticity  is  shed,  the  rate 
of  shedding  at  the  next  time  step  will  be  larger.  If  too  large  an 
amount  of  vorticity,  is  shed,  the  rate  of  shedding  will  be  smaller. 
This  kind  of  feedback  mechanism  suggests  the  existence  of  stability 
criteria  relating  the  grid  cell,  wing  panel,  and  time  step  size.  We 
have  not  yet  tested  the  possible  dependence  of  gross  predictions,  such 
as  aerodynamic  coefficients,  on  the  wing  panel izati on. 
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4.  Computational  Results  for  Two-dimensional,  Unsteady, 

Separated  Flow 

a.  Impulsively  Started  Flat  Plate  at  90  degrees  Angle  of  Attack 

The  first  problem  for  which  we  exercised  LIFTER  was  a  flat  plate  at 
90  degrees  angle  of  attack.  The  chosen  grid  was  40  by  40  cells,  and  the 
flat  plate  was  represented  by  22  panels  of  equal  length. 

Figure  4.1  presents  the  normal  force  coefficient  history  as 
predicted  by  LIFTER  compared  with  the  measured  average  value  (Hoerner, 
1975,  p.  21-1)  and  with  the  prediction  of  the  Helmholtz  model  (op. 
cit.).  There  are  no  time  dependent  measurements  of  the  force  history  on 
a  flat  plate  that  are  known  to  us.  Beyond  a  large  starting  transient  in 
the  normal  force,  there  is  no  evidence  of  oscillatory  behavior  in  the 
force  coefficient  that  would  be  indicative  of  the  formation  of  a  vortex 
street.  In  fact,  for  the  time  simulated,  the  force  coefficient 
decreases  monotonical ly  with  time,  as  would  be  indicative  of  the 
formation  of  a  stable  recirculation  cell. 

Figure  4.2  presents  a  series  of  comparisons  of  computed  streamlines 
with  aluminum  particle  streakline  photographs  and  Taneda  and  Honji  1971. 
We  can  expect  the  comparison  between  instantaneous  streamlines  and 
streak  lines  to  be  reasonable  if  the  exposure  time  of  the  photograph  is 
short  compared  to  a  characteristic  time  for  the  unsteadiness  of  the 
flow.  Taneda  and  Honji  do  not  state  the  exposure  time  for  their 
photographs,  although  we  may  presumably  infer  that  it  is  less  than  the 
spacing  of  their  photographs  in  time,  typically  0.5  B/U.  The  blockage 
factor,  that  is  the  width  of  the  plate  over  that  width  of  the  test 
channel,  was  less  than  10%. 

The  qualitative  behavior  of  the  computed  streamlines  compares  well 
with  the  streak  line  photographs  over  the  entire  series  of  comparions, 
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but  the  recirculation  cell  from  the  LIFTER  simulation  grows  signifi¬ 
cantly  more  slowly  than  the  recirculation  cell  in  the  streakline  photo¬ 
graphs.  Figure  4.3  makes  a  more  quantitative  comparison.  Taneda  and 
Honji  have  inferred  the  relation 

s/c  =  0.89(Ut/c)2/3 

where  s  is  the  length  of  the  recirculation  cell  and  c  is  the  chord 
of  the  plate.  Figure  4.3  compares  this  relation  with  the  length  of  the 
recirculation  cell  from  LIFTER  as  inferred  by  taking  twice  the  distance 
from  the  plate  of  the  centroid  of  the  vorticity  in  half  the  recircu¬ 
lation  cell.  The  prediction  from  LIFTER  does  not  follow  the  simple 
power  law  relationship  inferred  by  Taneda  and  Honji,  but  rather  grows  at 
first  more  quickly  and  then  more  slowly  than  the  inferred  experimental 
relationship.  Figure  4.3  also  contains  results  for  an  oscillating 
chord,  to  be  discussed  in  section  C. 

Both  the  photographs  and  the  computed  streamlines  confirm  the 
existence  of  a  stable  recirculation  cell  over  the  entire  length  of  the 
LIFTER  simulation.  For  large  Reynolds  numbers,  we  infer  from  figure  10 
of  Taneda  and  Honji  an  inviscid  scaling  for  the  critical  time: 

ter  =  5.3(c/U)  . 

While  the  determination  of  this  time  from  streakline  photographs  obvi¬ 
ously  involves  some  subjective  judgment,  we  may  expect  that  the  LIFTER 
computation  was  simply  too  short  to  exhibit  breakdown  of  the  recircu¬ 
lation  cell  and  consequent  rebounding  of  the  force  coefficient  from  its 
monotonic  decline.  It  also  is  not  clear,  from  the  paper  of  Taneda  and 
Honji,  whether  the  initial  breakup  of  the  recirculation  cell  is  a  two  or 
three  dimensional  phenomenon. 
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It  would  have,  of  course,  been  desirable  to  have  a  longer 
simulation,  but,  because  of  the  Kelvin  condition  and  because  of  the 
moment  calculation,  both  of  which  involve  integration  over  the  entire 
field  of  vorticity,  we  are  currently  forced  to  stop  our  simulations  when 
significant  vorticity  leaks  from  the  edge  of  the  computational  domain. 
We  will  discuss  the  economics  of  these  computations  further  in  Chapter 
Five. 


Figure  4.4  presents  plots  of  the  bound  vorticity  as  computed  by 
LIFTER  at  various  times.  Of  these,  the  only  one  that  is  readily 
interpretable  is  at  t  =  0,  where  the  exact  solution  is  known  for  the 
bound  vorticity  on  a  flat  plate  in  uniform  flow.  This  result  is  derived 
in  Appendix  B.  The  distribution  of  bound  vorticity  computed  by  LIFTER 
agrees  well  with  the  theoretical  prediction  for  all  but  the  last  two 
panels,  where  the  assumed  constant  distribution  of  vorticity  is  a  poor 
approximation  to  the  square  root  singularity.  It  can  also  be  seen  from 
figure  4.4  that  the  separated  vortex  flow  almost  inmediately  reduces  the 
bound  vorticity  to  a  small  fraction  of  its  initial  value.  The  residual 
concentration  of  bound  vorticity  at  the  end  of  the  wing  is  of  some 
interest.  We  will  speculate  as  to  the  meaning  of  this  phenomenon  in  the 
last  section  of  this  chapter. 

The  conclusions  we  may  draw  from  this  first  comparison  of  the 
predictions  of  LIFTER  with  experimental  data  are  not  as  strong  as  we 
would  like.  The  qualitative  behavior  of  the  code  predictions  is 
certainly  reasonable.  The  quantitative  behavior  of  the  recirculation 
cell  is  at  least  not  grossly  inconsistent  with  the  results  of 
measurements.  We  would  most  like  to  make  a  comparison  between  the 
measured  and  predicted  average  normal  force  coefficients,  but  this 
comparison  will  require  a  much  longer  computation  than  the  one  presented 
here. 
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b.  Impulsively  Started  Flat  Plate  at  45  Degrees  Angle  of  Attack 

Despite  the  fundamental  importance  of  the  flat  plate  at  90  degrees 
angle  of  attack  to  an  understanding  of  post-stall  aerodynamics,  the 
problem  is  peculiar  because  of  its  symmetry.  In  order  for  a  vortex 
street  to  appear,  it  is  necessary  for  some  symmetry-breaking  instability 
to  develop.  For  a  flat  plate  at  any  other  angle  of  attack,  asymmetric 
shedding  will  begin  immediately,  so  that  the  formation  of  a  long-lived 
recirculation  cell  is  unlikely. 

Figure  4.5  presents  a  sequence  of  vorticity  maps  from  our  next 
exercise  for  LIFTER,  a  flat  plate  at  45  degrees  angle  of  attack.  It 
seems  a  shame  to  belabor  the  elegance  of  the  simple  physics  in  this 
figure  with  technical  comment.  It  must  be  noted  that  the  aspect  ratio 
of  the  plots  is  not  quite  right;  the  dots  along  the  border  of  the  plots 
should  be  equally  spaced,  but  they  are  not. 

Figure  4.6  presents  another  striking  plot  from  the  same  simulation, 
this  time  of  the  absolute  value  of  vorticity,  and  this  time  with  the 
aspect  ratio  correctly  represented.  Next  to  it  is  an  aluminum  flake 
flow  visualization  of  the  same  problem,  reproduced  from  Cantwell  (1981). 
Although  the  photograph  is  obviously  in  a  much  later  state  of 
development,  the  similarity  of  flow  patterns  is  gratifying.  The 
shedding  frequency  in  the  numerical  simulation  is  obviously  very  close 
to  the  shedding  frequency  of  the  real  flow. 

Figure  4.7  presents  the  force  coefficient  history  from  this 
simulation,  and  here  we  begin  to  see  the  kind  of  effects  we  expect  from 
the  formation  of  a  vortex  street.  Comparison  of  the  force  coefficient 
history  with  the  flow  maps  shows  that  the  shedding  process  has  completed 
just  slightly  more  than  one  complete  cycle  at  the  end  of  the  simulation 
where  the  normal  force  coefficient  is  rising. 
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Figure  4.5.  Distribution  of  yorticity  for  an  impulsively  started 
flat  plate  at  45°  angle  of  attack. 
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ALUMINUM  FLAKE  FLOW  VISUALIZATION 
(Cantwell,  1981) 


Figure  4. 


Flat  plate  at  45°  angle  of  attack.  Comparison  of  computed 
distribution  of  absolute  value  of  vorticity  with  aluminum 
flake  flow  visualization. 
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Figure  4.7  also  shows  some  unexpected,  and  --  perhaps  -- 
disquieting  physics  in  the  form  of  a  significant  tangent  force!  The 
calculation  having  been  done  with  an  ostensibly  inviscid  code,  the 
circumspect  numerical  modeler  would  perhaps  have  done  better  to  have 
left  this  curve  off  the  plot,  but  there  it  is,  and  it  deserves  some 
explanation.  First,  of  course,  is  the  possibility  that  this  prediction 
is  simply  the  result  of  some  modeling  inaccuracy.  After  all,  despite 
the  intuitively  obvious  fact  that  the  tangent  force  on  a  flat  plate  in 
inviscid  flow  should  be  zero,  that  this  should  be  so  in  the  numerical 
simulation  requires  the  vanishing  of  a  time  derivative  of  a  particular 
moment  of  vorticity,  with  no  obvious  symmetry  in  the  problem  co  enforce 
this  condition. 

It  is  a  prediction  of  ideal  fluid  mechanics  that  there  is  no  net 
force  whatever  on  a  closed  body  in  a  uniform  freestream.  In  order  to 
explain  lift  and  drag  and  in  order  to  predict  them  with  LIFTER,  we  must 
introduce  viscous  mechanisms,  by  whatever  circuitous  route.  In 
classical  aerodynamics,  the  magical  device  is  the  Kutta  condition,  in 
LIFTER,  it  is  the  shedding  mechanism.  The  theorems  used  to  compute 
forces  in  LIFTER  remain  valid  if  viscosity  is  included,  and  LIFTER 
contains  mechanisms  that  are  manifestly  viscous.  There  is,  then,  no 
good  reason  to  dismiss  the  predicted  tangent  force  as  either  physically 
or  mathematically  incorrect. 

The  distribution  of  bound  vorticity,  presented  in  figure  4.8,  shows 
no  new  surprises.  The  vorticity  at  the  trailing  edge  bounces  up  and 
down  with  the  cyclic  shedding  process.  As  in  the  plate  at  90  degrees 
angle  of  attack,  there  remains  a  large  concentration  of  bound  vorticity 
at  the  leading  edge  throughout  the  simulation. 


Page  45 


PANEL  DIVISION 


igure  4.8.  Chordwise  distribution  of  bound  voriticty  for  an  impulsively 
started  olate  at  45°  anale  of  attack. 
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c.  Impulsively  Started  Flat  Plate  with  Oscillating  Cross  Section 
at  90  Degrees  Angle  of  Attack 

The  first  two  problems  in  this  section  are  important  from  the  point 
of  view  of  understanding  and  validating  the  basic  model.  The  underlying 
motivation  for  the  current  program  is  the  desire  to  achieve  sufficient 
active  control  of  the  shedding  process  as  to  have  significant  control 
over  the  net  forces  on  a  stalled  aerodynamic  shape.  We  felt  that  it  was 
desirable  to  get  results  as  early  as  possible  that  would  indicate  the 
feasibility  of  controlling  the  aerodynamics  through  manipulation  of  the 
salient  edges,  whether  through  oscillating  jet  flaps  or  oscillating 
mechanical  flaps. 

With  this  end  in  mind,  we  have  performed  two  simulations  of  flat 
panels  with  oscillating  chord  at  ninety  degrees  angle  of  attack.  The 
Strouhal  number  (the  ratio  of  the  oscillation  frequency  to  the  convec¬ 
tive  time  scale  based  on  chord)  was  chosen  to  be  two  for  one  of  the  runs 
and  0.5  for  the  other.  The  amplitude  of  oscillation  in  both  cases  was 
10%  of  the  chord  (that  is,  the  maximum  total  chord  is  1.1  times  the 
minimum  total  chord). 

The  results  of  these  simulations  are  presented  in  figures  4.9-4.10. 
The  normal  force  coefficient  history  for  the  Strouhal  number  2  case, 
presented  in  figure  4.10  is  in  striking  contrast  to  the  Strouhal  number 
0  case  previously  presented  in  figure  4.1.  Not  only  the  peak  normal 
force  coefficient,  but  also  the  mean  normal  force  coefficient  is  signif¬ 
icantly  enhanced  by  the  oscillating  chord.  The  normal  force  coefficient 
history  for  the  Strouhal  number  one  half  case  also  shows  augmented  mean 
lift,  but  the  augmentation  is  considerably  less  striking  than  for  the 
Strouhal  number  2  case. 
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Figure  4.10.  Normal  force  coefficient  history  for  an  impulsively  started 
with  oscillating  chord.  Strouhal  number  =  Z. 
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Referring  back  to  figure  4.3,  we  also  see  that  the  recirculation 
cell  for  the  plate  with  oscillating  chord  grows  faster  than  for  the 
plate  with  fixed  chord. 

We  would  expect,  even  for  an  unseparated  flow,  a  significant  effect 
from  oscillating  the  chord,  if  only  from  acceleration  and  deceleration 
of  the  changing  virtual  mass  associated  with  the  changing  chord.  The 
virtual  mass  associated  with  a  flat  plate  may  be  readily  inferred  from 
Batchelor  (1967  )  p.  431,  equation  (6.6.19)  to  be  ttpc^/^.  The  force 
associated  with  the  oscillation  of  the  virtual  mass  is  given  by 

3  1 

—  (irpc2u/4)  =  -  pUir2fc2(2A  sin(2irft) 

3t  2 

+  2A2  cos(2Trft)  sin(2TTft) 

The  normal  force  coefficient  associated  with  the  virtual  mass  is  then 
given,  for  small  amplitude  oscillations,  by 

2Tr2(fc/U)  A  sin(2irft)  =  2^2  Str  A  sin  (2TTft) 

which,  for  a  Strouhal  number  of  two  yields  an  envelope  amplitude  of 
roughly  4  with  an  amplitude  of  0.1,  and  an  envelope  amplitude  of  roughly 
2  for  the  same  amplitude  and  a  Strouhal  number  of  0.5.  These  numbers 
are  not  necessarily  valid  when  the  flow  is  separated,  as  in  our 
calculations,  but  they  agree  reasonably  well  with  the  envelope  of  the 
oscillating  component  of  the  computed  normal  force  coefficient. 

There  appears  to  be  significant  (greater  than  fifty  percent) 
augmentation  of  the  mean  normal  force  coefficients  for  the  cases  with 
the  oscillating  chord  as  compared  to  the  plate  with  fixed  chord.  The 
mean  for  all  three  cases  has  dropped  at  a  non-dimensional  time  of  6 
below  the  measured  average  normal  force  coefficient  for  a  flat  plate  at 
90  degrees  angle  of  attack. 
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d.  Conclusions  from  Two-Dimensional  Simulations 

Having  presented  our  results,  let  us  now  recapitulate  and  examine 
the  evidence.  Two  questions  should  be  answered: 

1.  Is  the  model  reasonable  in  the  sense  that  it  adequately 
predicts  reality? 

2.  What  have  we  learned  about  separated  flow? 

In  answering  the  first  question,  we  must  admit  that  the  evidence  at 
our  disposal  is  not  entirely  satisfactory  in  that  the  runs  are  shorter 
than  those  we  would  like  to  be  able  to  examine.  Only  in  one  case,  the 
flat  plate  at  45  degrees  angle  of  attack,  do  we  see  the  kind  of  periodic 
shedding  that  we  expect  is  responsible  for  the  failure  of  potential  flow 
and  Helmholtz  models  to  predict  the  measured  drag  coefficient  of  a  flat 
plate.  Even  in  that  case,  we  have  barely  one  full  cycle  of  the  periodic 
process  to  examine,  a  situation  that  makes  it  difficult  to  draw 
conclusions  about  the  mean  behavior  of  the  system.  Nevertheless,  it  is 
easy  to  believe  that  the  normal  force  coefficient  in  figure  4.7  is 
oscillating  about  the  correct  mean  value. 

The  qualitative  agreement  of  the  computed  wake  of  the  flat  plate  at 
forty  five  degrees  angle  of  attack  is  certainly  evidence  that  the  model 
not  only  predicts  the  periodic  shedding  process,  but  does  so  correctly. 

The  only  time  dependent  data  available  to  us  are  the  flow 
visualization  experiments  of  Toneda  and  Honji.  For  the  flat  plate  at 
normal  incidence,  the  persistence  of  a  stable  recirculation  cell  is 
evident  in  the  flow  visualizations  of  Toneda  and  Honji,  but,  based  on 
inviscid  scaling,  the  computed  recirculation  cell  grows  more  slowly 
longer  than  in  the  experiments. 
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Some  features  of  the  predictions  are  not  necessarily  contradicted 
by  observation,  but  certainly  run  contrary  to  expectations.  The 
prediction  of  significant  tangent  force  coefficient  for  the  flat  plate 
at  forty  five  degrees  angle  of  attack  must  be  included  in  this  category, 
as  must  be  the  not  unrelated  phenomenon  of  a  persistent  leading  edge 
singularity  for  this  calculation. 

We  have  learned  from  this  series  of  simulations  that  a  model  of 
vortex  shedding  that  ignores  the  details  of  the  boundary  layer  is 
capable  of  making  reasonable  predictions  of  the  shedding  process. 
Unlike  some  other  models,  our  model  does  not  predict  instantaneous 
relief  of  edge  singularities,  but  rather  that  reversed  flow  persists  for 
some  time. 

We  have  also  learned  that  periodic  forcing,  such  as  might  bo 
induced  by  jet  flaps,  enhances  both  the  rate  of  shedding  and  the  mean 
lift  coefficient. 
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5.  Directions  for  Further  Work  in  Two  Dimensions 
a.  Code  Modifications  Aimed  at  Economics 

We  have  been  limited,  in  our  first  series  of  simulations,  by  the 
need  to  capture  the  entire  field  of  shed  vorticity.  This  need  has 
arisen  both  because  of  the  Kelvin  condition,  and  because  of  our  use  of 
the  vorticity  theorems  of  Chapter  1  for  the  computation  of  aerodynamic 
forces.  We  can  easily  keep  track  of  vorticity  leaving  the  computational 
domain  for  purposes  of  enforcing  the  Kelvin  condition.  We  cannot  keep 
track  of  the  x  and  y  moments  of  vorticity  leaving  the  computational 
domain  without  further  assumptions,  such  as  assuming  that  the  vorticity 
simply  convects  with  the  free  stream  velocity.  Such  modifications  will 
be  worthv/hile  only  if  we  wish  to  make  very  long  simulations. 

The  cost  of  running  the  LIFTER  is  currently  dominated  by  the  cost 
of  solving  the  linear  system  of  equations  for  bound  vorticity.  The 
solver  currently  in  use  for  this  purpose  uses  Gauss  elimination  with 
pivoting.  Simply  switching  to  a  more  highly  optimized  solver  will  be 
the  most  cost  effective  step  we  can  make  in  obtaining  longer  urns  at  an 
acceptable  cost. 

b.  Calculations  Aimed  at  Further  Understanding  of  the 
Computational  Model 

We  have  so  far  been  unable,  because  of  costs,  to  check  the 
dependency  of  results  on  wing  panel i zation.  We  shall  surely  wish  to 
perform  such  an  experiment  when  more  efficient  software  for  computing 
bound  vorticity  is  in  place. 

We  may  also  wish  to  check  the  dependency  of  results  on  the  details 
of  the  vortex  shedding  model  by  ignoring  voriticty  swept  off  the  lee 
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side  when  the  flow  at  the  edge  is  outward  on  both  the  windward  and 
leeward  side  of  the  wing. 

c.  Calculations  Aimed  at  Further  Understanding  of  the  Physics 

and  of  the  Model 

We  would  of  course,  like  to  see  the  model  duplicate  average  drag 
coefficients  for  flat  plates,  and  to  do  this,  we  expect  to  need  much 
longer  runs  than  those  presented  in  Chapter  4.  Toneda  and  Honji  mention 
that  the  breakdown  of  the  recirculation  cell  is  an  unsymmetri cal 
process.  Such  an  asymmetry  must  correspond  to  an  instability  of  some 
kind.  Whether,  when,  and  how  such  an  instability  might  arise  in  a  two- 
dimensional  numerical  experiment  present  interesting  question.  The  free 
vortex  sheet  may  also  exhibit  a  symmetrical  instability. 

We  shall,  of  course,  want  to  simulate  a  slender  delta  wing  at  angle 
of  attack.  These  simulations  will  provide  us  with  a  further  opportunity 
to  compare  the  model  with  measurements. 

The  ultimate  goal,  of  course,  is  to  obtain  predictions  that  are 
relevant  to  active  vortex  control.  We  have  made  tentative  steps  in  this 
direction  with  the  runs  reported  in  Chapter  4,  but  we  shall  surely  want 
more  extensive  results  in  this  area. 
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6.  Numerical  Methods  for  Three  Dimensional  Unsteady  Flow 

The  two  dimensional  work  presented  so  far  is,  more  than  anything 
else,  preliminary  work  aimed  at  implementing  a  three  dimensional  model. 
The  understanding  we  have  achieved  as  a  result  of  the  two  dimensional 
work  reported  here  has  carried  us  a  considerable  distance  in  that 
di recti  on . 

The  ingredients  for  a  three  dimensional  model  are  much  the  same  as 
the  ingredients  for  a  two  dimensional  model;  an  Euler  solver,  a  means 
for  enforcing  zero  normal  velocity  on  embedded  surfaces,  and  a  means  for 
shedding  vorticity. 

Euler  solvers  for  incompressible,  three-dimensional  flow  are 
considerably  less  common  than  for  incompressible,  two-dimensional  flow. 
Poseidon  Research  has  recently  developed  one  such  code,  and  thic  code 
appears  to  be  well  suited  as  a  starting  point  for  a  three-dimensional 
version  of  LIFTER. 

Like  LIFTER,  Burner  uses  flux-corrected  transport  for  the  advection 
step.  Unlike  LIFTER,  Burner  is  formulated  with  three  components  of 
velocity  as  the  fundamental  field  quantities.  The  fact  that  Burner  does 
not  explicitly  conserve  vorticity  works  to  our  advantage.  In  partic¬ 
ular,  flows  simulated  by  Burner  will  separate  without  the  inclusion  of 
an  explicit  shedding  mechanism.  Figure  6.1  presents  a  field  of  flow 
vectors  obtained  by  simulating  a  square  jet  in  a  square  channel  with  a 
square  obstruction  midway  down  the  channel.  The  square  obstruction  is 
simulated  by  adding  impulsive  forces  to  the  model  system.  The  flow 
separates  both  at  the  inlet  and  at  the  square  obstruction. 

We  would  simulate  the  presence  of  eirtbedded  surfaces  for  aerodynamic 
calculations  in  much  the  same  way.  We  would  calculate  impulsive  forces 
that  cancel  the  normal  velocity  at  the  embedded  surface.  Just  as  for 
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bound  vorticity,  the  effect  of  the  impulsive  forces  is  non-local,  and 
calculation  of  these  forces  requires  the  solution  of  a  set  of  linear 
equations.  Application  of  the  local  effect  of  the  impulses  will  produce 
a  velocity  field  that  is  neither  incompressible  nor  zero  at  the  embedded 
surface.  The  solution  is  completed  by  solving  a  Poisson  equation  for  a 
velocity  potential  whose  gradient,  when  added  to  the  velocity  field, 
will  result  in  a  velocity  field  with  all  of  the  desired  characteristics. 

It  is  a  pleasant  circumstance  that  most  of  this  technology  is 
already  in  place,  and  that  the  only  new  software  required  will  be 
software  to  solve  for  the  panel  impulses. 

Although  we  are  confident  that  embedded  surfaces  will  shed 
vorticity,  it  is  by  no  means  clear  that  they  will  do  so  correctly. 
Exploring  the  actual  mechanism  by  which  the  code  produces  shed 
vorticity,  and  its  correctness,  will  constitute  a  significant  part  of 
the  development  effort. 
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Appendix  A 

Derivations  of  Force  and  Moment  Theorems 

The  two-dimensional  force  and  moment  equations,  (2.14)  and  (2.15), 
are  simple  enough  to  derive  in  Cartesian  component  notation.  The 
corresponding  derivations  in  three  dimensions  require  the  use  of  index 
notati on . 


a.  Two  Dimensions 


In  two  dimensions,  the  vorticity  equation  (2.11)  becomes 


3ai 

agy 

a9x 

—  +  U  *  Vo) 

s  - 

“  “  > 

at 

9x 

ay 

since  the  "vortex  stretching"  term  o)  •  Vu  is  zero  because  of  the 
absence  of  gradients  in  the  z-direction. 


Multiply  by  x  and  integrate  over  x  and  y: 


—  /  xoj  dA  +  /  (xu —  +  XV — 

dt  J  J  \  dx  ay 


j  dA 

ay  ' 


The  r.h.s.  can  be  rewritten  and  integrated  to  the  (distant)  boundary 
of  A: 
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3 

gy  —  (><gx) 
ay 


dA 


+ 


xg^Hy )  dC 


'>Jhere  n  is  the  unit  outward  normal  from  the  boundary  C.  The  second 
term  on  the  left  can  be  reworked  in  a  similar  fashion: 


3u  3  3v  j 

(xuoj)  -  xai - uo)  -  —  (xvoj)  -  X(jj —  dA 

3x  3y  Sy  1 


=  -  J"  m  dA  +  J"  (xutonx  -  xvwriy)  dC  , 


where  use  has  been  made  of  the  continuity  equation  3u/3x  +  3v/3y  =  0. 


The  perimeter  C  can  be  made  large  enough  so  that  g  and  u  are 
both  zero  on  it  {  oi  is  only  avected  by  the  fluid;  it  does  not  propagate 
instantaneously  like  a  potential  or  even  fast  like  a  sound  wave).  The 
line  integrals  are  then  zero,  and  we  have 
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But  0)  =  5V/3X  -  3u/3y,  SO  the  second  term  on  the  r.h.s.  becomes 


9u  I 

u — J  dA 

ay) 


3u  ) 

u—  dA 

ay ) 


Since  u  and  v  fall  to  zero  at  least  as  fact  as  r-1,  the  contour 
integral  is  again  zero. 

We  finally  have 


from  which 


xo)  dA 


A 


d  r 

P —  /  xu)  dA  . 

dt  J 

A 
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This  proves  the  y-component  of  (2.14).  Proof  of  the  x-component 
involves  multiplying  (Al)  by  y  and  proceeding  as  before. 

The  definition  of  moment  per  unit  leng'.h  in  two-dimensional  flow 
follows  from  (2.6): 


Me^ 


yfx)ez  dA 


(A2) 


To  prove  the  moment  relation  (2.15),  we  multiply  (Al)  by  r^  =  x^  +  y2 
and  proceed  as  before: 


dr  r  I  dui  V 

—  /  r2(jj  dA  +  /  (r^u —  +  r^v — )  dA 

dt  J  '  ax  9y' 

A  A 


=  f  .  r2~)dA  .  (A3) 

y  '  ax  dy 
A 


Dropping  contour  integrads  as  they  come,  we  find  for  the  r.h.s.  that 


a  1 

2xgy - (r2gx)  -  Zyg^  dA 

ay  ) 


(A4) 
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The  second  term  on  the  left  becomes 


—  (r^uoi)  +  — (r^voj)  -  r2a)j —  +  — )_  2xua)  -  2y 
3x  3y  \3x  3y/ 

A 

=  -  2  J"  (xm  +  yvoj)  dA 


vo)  dA 


/(  3v  3u  3v  3u) 

XU - XU —  +  yv—  -  yv— 

(  3x  3y  3x  3v ) 


3y  3x  3y  i 


/(  9  3u  9  /  u^\ 

j  — (xuv)  -  UV  -  XV - (x-  ) 

(Sx  av'  0  ' 


w  »  a  /  V 

- (yuv)  +  UV  -  yu —  +  — (y- 

3v  3y  3x'  2 


3x  3y'  2 


/(  9v  3u 

,xv —  +  yu  — 

'  9y  9x 


/(  3  v2  9  ,  u2 J 


9y  2  ax'  2 
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We  have  repeatedly  made  use  of  the  assumption  that  quanities  like  xv2, 
yu^,  etc.  fall  to  zero  factor  than  r"^.  That  is  so  by  virtue  of  (2.20), 
which  implies  that  velocity  falls  as  r"2  far  from  the  region  occupied 
by  vorticity. 

Combining  (A2)-{A5),  we  obtain 

P  d 

M  =  -  — 

2  dt 

which  is  the  same  as  (2.15). 


b.  Three  Dimensions 

The  proofs  of  the  two-dimensional  force  and  moment  theorems  are 
fairly  complicated.  Matters  get  out  of  hand  in  three  dimensions  unless 
index  notation  is  used.  With  index  notation,  the  i^  components  of 
velocity  and  vorticity  are  represented  as  ui  and  a)i  respectively. 
The  curl  involves  the  permutation  symbol  Eijk>  so  that 

auk 

ij^i  “  ^ijk -  • 

9Xj 

eij|(  is  +1  if  ijk  are  an  even  permutation  of  123,  -1  if  ijk  are 

an  odd  permutation  of  123,  and  zero  if  two  of  the  indices  are  the  same. 
Another  symbol  of  great  utility  is  the  Kronecker  delta  djj,  which  is 
+1  if  i  and  j  are  equal  and  zero  if  not.  The  continuity  equation 
can  thus  be  written  as 


axj 


(A7) 
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Whenever  indices  are  repeated,  as  in  (A6)  and  (A7),  they  are  summed  over 
1,  2,  3  (the  Einstein  sunmation  convention). 

e  and  6  are  connected  by  the  identity 

£ijk  £imn  “  Sjm  Skn  “  ^km  5jn  »  (A8) 

used  frequently  in  the  following  proofs.  Often  we  must  rearrange  the 
indices  of  e  to  apply  (A8).  From  the  definition  of  e, 

^ikj  “  '  £ijk  »  ekij  “  ^ijk  •  (A9) 

Odd  and  even  permatations  of  e’s  indices,  exemplified  by  (A6),  are 
used  in  the  following  without  comment. 

In  index  notation,  the  three-dimensional  force  theorem  (2.14) 
becomes 


P  d 
2  dt 


xj  wk  dV 


(AlO) 


where  the  first  equality  follows  from  the  definition  (2.5)  of  ^  (with 
^  =  p£),  and  the  second  is  to  be  proved. 


The  proof  is  based  on  the  vorticity  equation  (2.11),  which  becomes 


9(Dk  9ajk  BUk  Bgp 

—  +  um - u)m —  =  Ekmn — 

3t  BXni  BX|n  BX|„ 
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in  index  notation.  Take  the  cross-product  with  x: 


3U|(  \  9gn 

^ijk  ^ijk  (“ffl  “  wm  )  =  ^ijk  ^kmn  '  (All) 

3t  9Xjn  3x^1 '  Sxp, 


Apply  the  e-6  identity  to  the  r.h.s.: 


99n 

39n 

^ki j  ^kmn 

"  (*Sim  5jn  ■  'Sjn  <5in)xj 

9Xm 

9xm 

9gj  9gi 

=  - ^3  — 

9x-| 

9  9 

“  -  (^jOj)  -  9j6ji  -  —  (><j9i)  9ii5jj 

9Xi  9xj 


9  9 

=  —  (xj-gi) - (xjgi)  +  Zg^  ,  (A12) 

9x-j  9xj 


where  we  have  used  the  facts  that  9xj/9x-j  =  6  j  i  and  6jj 
1  +  1  +  =  3.  Rework  the  second  term  on  the  l.h.s.  in  like  fashion: 
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3wk  3Uk 

^ijk  (‘^n - t*Jni — 

’  9Xn)  3Xm 


(  3 

“  ^ijk  ‘^k  <Sjin 

(  3x^1 


-  (Xj(i)n,Uk)  -  u)(n  Uk  6jnn 

9Xff, 


■  —  pijk  "  ^^ijk  ‘^k  ’  (^13) 

3Xn)  (  ) 


Moreover , 

3un 

^ijk  <*>k  “  ^ijk  ekmn  — 

3Xn, 


=  ('SiiD'Sjn  ■  6jni6in)  — 

3Xm 


3Uj  3U'|  9  ^  j  ^ 

=  uj -  -  Uj — ■  =  —  -  -  —  (u-jUj)  {A14) 

3Xt  3Xj  3Xt  2  3Xj 

Gathering  (All)-(Ai4),  integrating  over  an  arbitrarily  large  volume 
V,  and  eliminating  surface  integrals  involving  ui  or  g^-  ,  we  obtain 
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Velocity  components  ui  fall  at  least  as  fast  as  r-2  in  three 
dimensions,  so  the  surface  integral  vanishes,  and  we  have  proved  (AlO). 


The  definition  of  moment  (2.6)  and  the  moment  theorem  (2.17) 
require  us  to  prove  that 


eijk  xj  gk  dV 


V 


1  d 
3  dt 


^ijk  Xj  (ekmnxmwn)  dV 


V 


(A15) 


The  proof  is  complicated  and  will  only  be  sketched  here.  Cross  the 
vorticity  equation  twice  with  x: 


^  n 

^ijk  ^kmn  %  —  + 

at 


c-  c  / 

^ijk  ^kmn  ^m  H^p  — 

^  axp 


CO  p 


axp 


9gq 

^ijk  ekmn  ^npq  -  •  {A16) 

axp 
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Manipulate  the  r.h.s.: 

I  39q  \ 

^ijk  l^nkm  ^npq  - / 

9Xp 


agq 

■  ^ijk  (<Skp  6fnq  ■  <Smp  6kq) 

3Xp 


/  agm  3gq 

■  ^ijk  % - — 


(  ^ 

"  ^ijk  —  {xjXnjgm)  -  Xj  6kni  •  gm 

(  axk 


^  I 

- Xj  9k  Snm  gk  iSnij  ( 

3Xk  > 


'  ^ijk  ) -  (xjXnigin)  -  -  (xjXg,gk) 

( 3Xk  axn, 


-  Xj  gk  +  3xj  gk  +  Xj  gkj 

The  second  term  on  the  left  of  (A16)  can  be  written  entirely  in  terms  of 
divergences.  Integrating  (A15)  over  all  space,  we  find  that 
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d 

dt 


J 

V 


Eknm  % 


=  3 

V 

where  the  surface  integrals  vanish  provided  that  ui-j  and  g,  are 
confined  to  a  finite  domain,  and  that  XjUmUn  falls  to  zero  faster  than 
r~2j  which  is  always  true.  This  confirms  (A15)  and  hence  (2.17). 
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Y(9  )  s1n9'  =  Aq  +  B  cose' 

«o  t,e  zero  to  satisfy  t«e  circulation  theorm.  Then 

u  s1n  Ot  =  Bo/2  . 


'hus 


1  -  2x/C 
V(x/cj[l-{x/c)] 


Y 


2  U  sinactn  9 


U  sin 


